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0^ Abstract. We present a method for finding the eigenmodes of the Laplace operator 

■ acting on any compact manifold. The procedure can be used to simulate cosmic 

microwave background fluctuations in multi-connected cosmological models. Other 
applications include studies of chaotic mixing and quantum chaos. 

1. Introduction 

Much of physics boils down to solving differential equations subject to certain 
boundary conditions. Waves in a box are a classic example. The box we have in 
' mind here is some closed manifold, and the waves are those of a massless scalar 

field. This picture arises naturally when studying density perturbations in a multi- 
connected universe^]. The mathematical problem can be stated: Find all square 
integrable functions \&g(x) that satisfy the partial differential equation 

(V + 9 2 )v]/ 9 (x) = 0. (1) 
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Here V is the Laplace operator on some closed manifold X and the constant q is 
an eigenvalue of the Laplacian. The complexity of the problem is controlled by 
the geometry of X. When S is n-dimensional Euclidean space, E n , modulo some 
discrete group of covering transformations, T, it is a simple matter to write down 
analytic expressions for the eigenmodes ^ q (x). In contrast, if S describes some 
compact hyperbolic manifold, H n /T, then the eigenmodes cannot be expressed in 
closed analytic form. This difficulty is closely related to the chaotic behaviour of 
geodesic flows on compact negatively curved spaces [^|. Here we describe a numerical 
solution to the problem based on Fourier filtering solutions to the scalar wave equation 
□ \E' = 0. While our approach works for any topology, we will focus on hyperbolic 
manifolds as these are of the most interest to cosmology. 

2. Ringing out the modes 

Rather than attack Laplace's equation directly, we begin by introducing a ficticious 
time dimension t, thereby lifting the problem to solving the wave equation 

on the Lorentzian manifold A4 = R x S. Consistent with this picture we adopt a 
product ansatz for the eigenmodes: 

#q(t, x) = exp(-iw g t)# g (x) . (3) 



2 



The eigenfrequency u q is fixed by (|^) to equal q. Assuming for now that the spatial 
boundary conditions have been properly enforced, we can evolve the initial data 
^(OjX), 9t^(t, x)|q = according to (|2|) to find \t(£,x). By Fourier transforming 
in time: 

/>oo 

a w (x)= / *(t,x)e Mt di, (4) 
Jo 

and calculating the power spectrum 

P(w)= / KlVstfa;, (5) 



s 

we are able to isolate the eigenfrequencies w g . These are located at local maxima of 
P(u>). Once the eigenfrequencies are known the individual spatial eigenmodes can be 
extracted: 

%(*)= lim I [ T ^(t,x)e^dt. (6) 

In practice the integration time T will be finite, so we can only resolve modes separated 
in frequency by at least Au = 2ir/T. 

The main difficulty in implementing the above procedure stems from the 
complicated periodic boundary conditions that are imposed by the topology. The 
remainder of this paper is devoted to describing a numerical solution to this problem, 
and illustrating how it works by finding the eigenmodes for a genus 2 surface with 
constant negative curvature. 

3. The double doughnut 

The manifold we will focus on is the two-hole doughnut S — H 2 /T, where T is a 
discrete subgroup of 50(2, 1) with presentation 

r = {So, gi, 52, 53 : 9ogi 1 9293 1 9o 1 9i92 1 93} ■ ( 7 ) 
The group generators have the SO(2, 1) matrix representation|| 

cosh rj sinh r\ 

cjj = Rj I sinh 77 cosh 77 | R~ x , (8) 
1 

where ij = 2 arccosh(l + a/2) and 

/ 1 
Rj=\ cos(i7r/4) - sin(j7r/4) ) . (9) 
y sin(j7r/4) cos(j7r/4) 

The double doughnut can be obtained by gluing together opposite faces of a regular 
octagon with dihedral angles 7r/4. The octagonal fundamental domain (FD) is drawn 
in Fig. 1 using the Poincare disc model for H 2 . 
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Figure 1. The fundamental domain for the double doughnut (shown bordered by a 
solid line) in the Poincarc disk (bordered by a dashed line). 

The double doughnut has area 4tt and diameter d = 2 arctanh(2~ 1 / 4 ) w 2.448. All 
distances arc quoted in units of the curvature radius. The largest simply connected 
circle that can be drawn inside S has radius r\- = arccosh(l + \/~2), and the smallest 
circle to fully enclose the fundamental domain has radius r] + = 2arccosh(\/ 2 + \/2). 
The light crossing time varies in the interval 2ry_ < t c < 2rj + , and has the mean value 
t\ = \Z4~7r w 3.545. This is the Lyapunov time for chaotic flows on S, and represents 
the characteristic dynamical timescale for our system. 



4. Solving the wave equation 

The dynamics is described by a massless scalar field </> with action 

s = J V^gg^d^d^d 3 x, 



(10) 



evolving on the Lorentzian manifold (R x S, g). Choosing a Poincare metric for H 2 
we have 



ds 2 



-dt 2 + 



Ajdx 2 + dy 2 ) 
(1 — x 2 — y 



2\2 " 



(11) 



The coordinate distance r = (x 2 + y 2 ) 1 ^ 2 is related to the proper distance rj by 
r = tanh(r//2). To evolve the system numerically we begin by discretising time and 
space such that x = i Ax, y = j Ay and t — k At. The discretised action then reads 



{<f>(i + l,j,k)-cf>(i,j,k)) 2 (0(i,j + l,fc)-0vU»*)) 



Ax 2 



+ 



Ay 2 

(cj)(i,],k + l)-cj)(i,j,k)) 2 



(1 - (iAxf - (jAy) 2 ) 2 



At 2 



Ax Ay At. (12) 
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The equations of motion that follow by varying with respect to </>(i, j, k) are: 

ASU(i,j, k) = (1 - (zAx) 2 - (jAy) 2 ) 2 [6 2 <j>(i,j, k) + 6 2 <j>(i,j, k)] (13) 



where 




(a) (b) 

Figure 2. (a) The grid points are first sorted into masters (dark crosses), slaves 
(heavy dots) and freemen (light crosses), (b) The slave images (heavy dots) are 
located inside the FD. 



Our next task is to enforce periodic boundary conditions at the edge of the 
fundamental domain. To do this we begin by sorting all points in the spatial grid into 
those inside and outside the fundamental domain. There is a simple sorting algorithm 
that works for all topologies. In what follows, |x| denotes the proper distance from 
the origin. 

Algorithm 4.0. Point Sorting 

(i) If |x| < ?7_ then the point lies inside the FD. If |x| > 77+ then the point lies 
outside the FD. 

(ii) For ry_ < |x| < ry + , act on x by all n group generators and their inverses to form 
the 2n image points x' ± i = gf 1 *-- 

(iii) If |x' ±i | < |x| for any i, then x lies outside the fundamental domain, else x is 
inside the fundamental domain. 

Only the "master" points inside the FD need to be evolved. However, forming the 
second derivatives S 2 <fi and 5 2 4> for the masters often requires knowledge of field values 
outside the FD. We refer to all points outside the FD that are within one grid point 
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of a master as "slaves". Points that are neither masters nor slaves play no part in 
the numerical evolution and arc designated "freemen" . Each slave has a unique image 
inside the FD. We can find the position of these images using the following algorithm. 

Algorithm 4.1. Locating the fundamental image 

(i) Act on x by all n group generators and their inverses to form the 2n image points 

(ii) Find the image point x' ± z nearest the origin and call it x'. 

(iii) If |x'| < |x| then let x = x' and go to (i), else x is the fundamental image. 

Since typical coordinate grids are not mapped into themselves by the fundamental 
group, the slave image points will not lie on the computational mesh. Therefore we 
have to interpolate to find the field value at each slave point. The procedure of point 
sorting and image finding is illustrated in Fig. 2. Only one quadrant of the FD is 
shown since the FD has 8-fold symmetry. The slave images are those of slave points 
located in the other three quadrants. 
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Figure 3. Interpolating the image of a slave point. 



The interpolation of the slave values is done using either a 3-point, 4-point or 
6-point interpolation scheme. Employing the labelling convention shown in Fig. 3, the 
6-point interpolation is described by 

4>{p Ax, q Ay) = q(q - l)/2 <fo-i + p(p - l)/2 <£_i 

+ (1 - pq - p 2 - q 2 )4> m + p(p -2q+ l)/2 <j> W 

+ q(q -2p+ l)/2 (j>oi + pq<t>u + 0(e 3 <M • (15) 

Here e = max(Ai, Ay). By rotating the coordinates so that the interpolation is 
performed in the quadrant x > 0, y > 0, we can ensure that the (0, 0) vertex of the 
cell enclosing the image point lies nearest the origin. This increases the chances that 
(—1,0) and (0,-1) are master points. Nevertheless, we need to check that all points 
involved in the interpolation are either masters or slaves. If a freeman is used in the 
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interpolation it must be enslaved. Once this is done for all the slave points we arrive 
at a set of N coupled linear equations involving N slaves and M masters (the slave 
drivers). Writing the list of slaves and slave drivers as the column vectors <j>s and 
4>SD, the interpolation equation (15) can be used to form the linear system 



A<j> s = 



= A B<j>s 



D ■ 



(16) 



where A is a N x N matrix and B is a N x M matrix. This procedure only has to be 
performed once at the beginning of a simulation. The N x M matrix C = A _1 B can 
then be stored and called during the numerical evolution. In summary, at each time 
step the masters are evolved according to (13). The slaves are then updated using 
(Pf). The slaves are the glue that holds together identified sides of the fundamental 
cell. If the glue is poor energy can either leak in or out of the system. Thus, by 
keeping track of the total energy E = K + G we can check that the periodic boundary 
conditions are being properly implemented. The kinetic energy is given by 



K k = 



E 



4(0(i,j, k)-<j>{i,j, k-l)f 



(1 - {iAx) 2 - {jAy) 2 ) 2 At 2 



Ax Ay . 



and the gradient energy is given by 



G k = £ 

i,3 



{4>{i + l,3,k)-(j){i,j,k)f , {(j>{i,j + l,k)-(j){i,j,k)) 



2 1 



Ax 2 



Ay 2 



(17) 



AxAy, (18) 



Results 



Several simulations were run using 301 2 grids covering the region —0.8 < x,y < 0.8. 
The coordinate grid spacings were: Ax = Ay = 2/375 and At = Ax/10. The 
proper distance between grid points varies between Sx — Ax near the origin and 
Sx = (2 + y/2) Ax at the extremities of the FD. Consequently, the resolution is ~ 3.4 
times worse at the edge of the FD. Other metrics could be used to give a more even 
coverage if desired. As a check on the accuracy of the code the area was evaluated 
and compare to the continuum answer: 



.4 



masters 

E ■ 



AAxAy 



1,3 



(1 - (iAx) 2 - (jAy) 2 ) 2 



(1.00014) 4tt. 



(19) 



Two types of initial conditions were used for 0(0, x). The first employed a random 
superposition of Gaussian peaks and was designed to excite as many eigenmodes 
as possible. The second employed a superposition of two spherically symmetric 
eigenmodes (q = 0.1 and q = 0.5) of the Laplacian on H 2 . This choice was designed 
to excite low lying eigenmodes. The waveforms were smoothed to remove any large 
gradients caused by the imposition of the periodic boundary conditions. The residual 
monopole contribution was removed to increase our chances of resolving any low lying 
multipoles. A stroboscopic sequence showing the evolution of the random Gaussian 
peaks is shown in Fig. 4. 

The wave was evolved for 2 17 timesteps, which equates to T = 69.9051 curvature 
units, or roughly 20 light crossing times. The total energy remained roughly constant 
throughout the simulation, although there were ~ ±2% fluctuations about the mean 
value. We believe these variations are due to the uneven grid resolution. 
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Figure 4. The evolution of an initially static waveform shown every 0.2 time units. 
Time runs cartoon style. 

The wave was Fourier analysed to produce the scaled power spectrum shown in 
Fig. 5. The power spectrum is multiplied by y/U for lo > to reflect the larger 
statistical weight of the high frequency modes. 
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Figure 5. The scaled power spectrum for the wave shown in Fig. 4. 

The first three peaks are at u)\ = qi = 1.97 ± 0.09, q2 — 2.30 ± 0.09 and q% = 
2.90 ± 0.09. The errors reflect the finite frequency resolution of Acj = 2-k/T — 0.09. 
Aside from some residual monopole power at q — there is no evidence for any low 
lying eigenmodes. We will return to this crucially important point later. 




Figure 6. The two lowest eigenmodes found by Fourier filtering. 

By Fourier filtering at u) = lu± and 10 = loi we are able to extract the corresponding 
spatial eigenmodes. These are displayed in Fig. 6. By numerically evaluating the 
Laplacian of these waveforms we were able to confirm that they are indeed eigenmodes 
with the correct eigenvalues. We illustrate this in Fig. 7 by displaying the eigenmodc 
<f> q3 (x.) and — q^ 2 V(f> q3 (x). The agreement is remarkably good over most of the FD, 
save for some small regions where the spatial gradients are large. The eigenmodes can 
be improved by a second scrubbing run. To do this we use the approximate eigenmode 
as initial data and evolve the system as before while filtering at the appropriate 
eigenfrequency. This helps to remove contamination from the other eigenmodes. 
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Figure 7. (a) The 53 eigenmode. (b) The q-j eigenmode acted on by — q 3 2 V. 

One possible explanation for the lack of eigenmodes below q\ = 1.97 might be a 
lack of long range power in our choice of initial conditions. To check this we tried 
a different initial waveform based on a superposition of the spherically symmetric 
q = 0.1 and q = 0.5 eigenmodes of H 2 . When the periodic boundary conditions are 
imposed the resulting initial waveform has 8-fold symmetry. 




Figure 8. The 8-fold symmetric q = 0.1,0.5 initial waveform and the resulting 
scaled power spectrum. 

The scaled power spectrum shown in Fig. 8 shows no sign of any modes below 
qi = 1.97. As expected, many of the modes excited by the random superposition 
of Gaussian peaks are missing from the spectrum produced by the 8-fold symmetric 
initial conditions. However, the modes that are present agree with those in seen in 
Fig. 5. 
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6. Supercurvature modes 

In n dimensions the eigenmode spectrum for H n takes all values in the range 
q = [<7o,oo), where qo = (n — l)/2. Modes with q < qo are not square integrable 
in infinite n-dimensional hyperbolic space, but they are square integrable in the 
compact quotients S™ = H n /T. According to BuserQ, many compact hyperbolic 
manifolds support supercurvature modes with q < q . The name supercurvature 
refers to the fact that these modes support significant long range correlations on scales 
larger than the curvature scale. In contrast, modes with q > qo have exponentially 
damped correlations outside the curvature radius, even though their wavelengths 
A = 2n(q 2 — f/g)^ 1 ^ 2 may greatly exceed the curvature scale. This suppression can 
be understood on the grounds of flux conservation in a space where the volume 
grows exponentially on large scales. Supercurvature modes are very important in 
a cosmological setting as even a single supercurvature mode can significantly increase 
the amplitude of cosmic microwave background fluctuations on large angular scales Q. 

There are a number of upper and lower bounds for the lowest eigenvalue, q\ > qo, 
of the Laplacian on H n /T . In two dimensions the tightest upper bound we know of 
was found by Cheng 0: 

*?<M^y< (20) 



4 V d , 

where d is the diameter of E. Applying this to our genus 2 example we find that 
qi < 2.61, which is consistent with what we found (qi = 1.97±0.09). Using Cheeger's 
inequality, q\ > h 2 /4, where h is Cheeger's isoperimetric constant, Balazs & Voros|| 
suggest that qi > 0.7. This is also consistent with our result. 

To be certain that we have not missed any low lying eigenvalues we need to 
use other methods such as the Selberg-Gutzwiller global eigenvalue count |7). After 
presenting our results in Cleveland, we discovered a paper by Aurich & SteinerQ 
that calculates all the low lying eigenvalues for the double doughnut using a suitably 
regularised Gutzwiller trace formula. They find the lowest three eigenvalues to be 
q% = 1.959, qi = 2.314 and qz = 2.872, in perfect agreement with our values. This 
proves that we did not miss any low lying eigenvalues. While their values are sharper, 
the advantage of our method is that it allows us to extract the eigenmodes in addition 
to the eigenvalues. 

The lack of low lying eigenmodes on the regular octagon is probably related to its 
high degree of symmetry. Indeed, it is possible to move to an alternative description 
of the manifold using a larger fundamental group and a fundamental cell 96 times 
smaller than the octagon^). Viewed from this perspective, the lowest eigenmode is 
remarkably low, as it corresponds to a wave with a wavelength ~ 10 times larger than 
the diameter of the desymmetrised fundamental cell. 



7. Concluding Remarks 

Having demonstrated that our approach works in 2-dimensions, our next task is 
to apply it to the cosmologically relevant case where E is a compact hyperbolic 
3-manifold. The computational cost of having an extra dimension will be offset 
by the availability of examples with very small fundamental domains. Many small 
hyperbolic 3-manifolds have FD's that have outradii, 77+, smaller than the radius 
of curvature ||. Consequently, space looks approximately Euclidean inside the FD, 
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so standard coordinate systems such as the Klein metric will produce fairly uniform 
computational grids. Once we have found all the low lying eigenmodes we can use 
them to produce simulated maps of the cosmic microwave background radiation. These 
maps can then be compared to observational data, or used to test proposals for finding 
the large scale topology of the universe px||. 
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